from scipy.interpolate import RegularGridInterpolator
from scipy.interpolate import interpn
import numpy as np

def value_func_3d(x, y, z):
    return 2 * x + 3 * y - z

if __name__ == '__main__':
    x = np.linspace(1,4,11)
    y = np.linspace(4,7,22)
    z = np.linspace(7,9,33)

    # x = np.linspace(0,5)
    # y = np.linspace(0,5)
    # z = np.linspace(0,5)

    pcs = (x, y, z)
    values = value_func_3d(*np.meshgrid(*pcs))

    data = value_func_3d( *np.meshgrid(x,y,z, indexing='ij', sparse=True))
    itp_f = RegularGridInterpolator((x, y, z), data)

    query_point = np.array([[2.1, 6.2, 8.3], [3.3, 5.2, 7.1]])
    # query_point = np.array([2.21, 3.12, 1.15])
    query_value_1 = itp_f(query_point)
    query_value_2 = interpn(pcs, data, query_point)

    print('Done!')